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ABSTRACT 

♦ 

In  this  paper,  we  will  present  a new  representation  of  a probability 
density  function  on  the  three  dimensional  rotation  group,  S0(3),  which 
generalizes  the  exponential  Fourier  densities  on  the  circle.  As  in  the 
circle  case,  this  class  of  densities  on  S0(3)  is  also  closed  under  the 
operation  of  taking  conditional  distributions.  Several  simple  multistage 
estimation  and  detection  models  will  be  considered  in  this  paper.  The 
closure  property  enables  us  to  determine  the  sequential  conditional  den- 
sities by  recursively  updating  a finite  and  fixed  number  of  coefficients. 
It  also  enables  us  to  express  the  likelihood  ratio  for  signal  detection 
explicitly  in  terms  of  the  conditional  densities. 

An  error  criterion,  which  is  compatible  with  a Riemannian  metric, 
will  be  introduced  and  discussed  in  this  paper.  The  optimal  orientation 
estimates  with  respect  to  this  error  criterion  will  be  derived  for  a given 
probability  distribution,  illustrating  how  the  updated  conditional  densi- 
ties can  be  used  to  sequentially  determine  the  optimal  estimates  on  S0(3). 


This  research  was  sponsored  by  the  Air  Force  Office  of  Scientific  Research, 
Air  Force  Systems  Command,  USAF,  Under  Grant  No.  AFOSR-74-2671B. 


I. 


Introduction 


Rigid  body  rotations  are  involved  in  many  important  practical  problems 
of  detection,  estimation,  and  control.  Some  notable  examples  can  be  found 
in  gyroscopic  analysis  and  satellite  attitude  determination  and  control. 

While  linearization  and  approximation  techniques  have  led  to  many  useful 
results,  simple  analytic  tools  which  will  enable  us  to  analyze  and  synthe- 
size the  optimal  structures  have  long  been  desired. 

Optimal  estimation  and  detection  schemes  for  discrete-time  processes 
whose  state  space  is  a circle  or  sphere  have  been  obtained  in  [1]  and  [2] 
by  using  a novel  representation  for  probability  densities  which  has  the 
form  exp  f where  f is  a finite  linear  combination  of  functions  which 
form  a complete  orthogonal  system  on  the  state  space  involved.  In  the  case 
of  the  circle,  circular  functions  were  used,  while  both  spherical  harmonics 
and  multiple  trigonometric  functions  were  employed  for  densities  defined  on 
the  sphere. 

In  this  paper  the  same  approach  will  be  taken  for  discrete-time  rotation- 
al processes  by  introducing  a similar  exponential  density  referred  to  as  a 
rotational  exponential  Fourier  density  (REFD)  defined  on  the  group  of  rota- 
tions of  three-dimensional  space,  that  is  obtained  by  using  a sequence  of 
irreducible  unitary  representations  which  form  a complete  orthogonal  system 
on  S0(3).  It  can  be  shown  that  a continuous  density  function  on  S0(3) 
can  be  approximated  by  such  a RF.FD  as  closely-  as  we  wish  in  the  soace 
of  square  integrable  functions. 

As  in  the  circle  case,  the  class  of  REFD's  of  a certain  finite  order  is 
closed  under  the  operation  of  taking  conditional  distributions  as  a consequence 
of  the  group  structure  of  S0(3).  We  note  that  there  does  not  exist  such  a 


closure  property  in  the  sphere  case  [2]  and  a combined  usage  of  two  kinds  of 
exponential  densities  is  required  to  treat  analogous  estimation  and  detection 
problems  on  the  sphere.  It  ^rill  become  clear  in  the  sequel  that  it  is  exactly 
this  closure  property  of  REFD's  that  yields  simple  estimation  and  detection 
schemes  which  update  the  sequential  conditional  densities  by  recursively 
revising  a finite  and  fixed  number  of  parameters. 

One  drawback  of  REFD's  is  that  there  is  no  known  closure  property  of 
convolution,  which  places  a restriction  on  the  type  of  rotational  signal 
processes  that  can  be  considered  for  the  above  approach.  More  specifically, 
the  rotational  signal  process  considered  in  this  paper  disallows  random  driving 
terms.  Given  two  orientations  of  a rigid  body  which  are  represented  by  two 
points  on  S0(3),  the  minimal  angle  in  radians  required  to  bring  one  into  the 
other  is  a Riemannian  metric  on  SO (3)  and  is  a natural  measure  of  the  distance 
between  them. 

An  error  criterion  for  orientation  estimation,  which  is  compatible  with 
this  measure  of  distance  will  be  introduced  in  this  paper.  It  is  compatible  in 
the  sense  that  it  is  an  increasing  function  of  this  Riemannian  metric.  Various 
descriptions  and  properties  of  this  criterion  will  be  discussed.  The  optimal 
orientation  estimate  and  its  estimation  error  with  respect  to  this  error 
criterion  will  be  derived  for  a given  probability  distribution,  thereby 
illustrating  how  the  updated  conditional  densities  can  be  used  to  sequentially 
determine  the  optimal  estimates  of  the  rigid  body  orientations. 

It  should  be  mentioned  that  estimation  for  continuous-time  directional 
processes  on  S0(3)  has  been  considered  in  [9]  and  [10]. 
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II.  Preliminaries  and  Rotational  Exponential  Fourier  Densities 


A rotation  of  three-dimensional  Euclidean  space  about  a fixed  point  can 
be  described  analytically  as  a linear  transformation  of  the  space  that 
preserves  distances,  leaving  the  origin  unchanged,  which  can  be  represented 
by  an  orthogonal  matrix.  In  this  paper  we  will  be  concerned  with  the  group 
of  rotations,  denoted  by  S0(3),  that  consists  of  those  rotations  whose  matrix 
representation  has  determinant  equal  to  +1;  that  is,  we  are  excluding  those 
rotations  that  consist  of  a rotation  followed  by  a reflection.  There  are 
several  different  ways  of  parametrizing  SO (3)  such  as  by  Euler  angles,  unit 
quaternions,  direction  cosines,  and  Cayley-Klein  parameters  which  are  all 
discussed  in  [3].  While  we  will  have  occasion  to  refer  to  all  of  these  methods, 
it  will  be  customary  to  parametrize  any  rotation  R belonging  to  SO (3)  by 
a triple  of  Euler  angles  (<J>,0,i{/)  which  determine  the  rotation  according  to 
the  following  sequence  of  rotations  [4,  p.l07J: 

(i)  a rotation  through  the  angle  , 0 £ 6 < 2tt  , about  the  z-axis, 

(ii)  a rotation  through  the  angle  0,0  ^ 0 < it  , about  the  new  x-axis, 

(iii)  a rotation  through  the  angle  ip,0  < ip  < 2ir  , about  the  new  z-axis. 

These  rotations  are  illustrated  in  Figure  1,  with  x,  y,  z being  the  final 

position  of  the  three  coordinate  axes. 


Fig.  1 
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An  orthogonal  matrix  R that  corresponds  to  this  rotation  can  be  obtained 
as  the  product  of  the  three  matrices  corresponding  to  the  rotations  (i) , (ii) , 
(iii):  R * Z(4>)X(0)Z(i^)  where 


cos  <f>  - sin  (f> 

7(<t>)  = sin  <fi  cos 


cos  0 
sin  0 


- sin  0 
cos  0 


so  that 

cos  (p  cos  \p 

- sin  <p  cos  0 sin  tp 

R = sin  <})  cos  ip 

+ cos  cj)  sin  ip  cos  0 

sin  sin  0 


- cos  <f>  sin 

- sin  <(>  cos  ip  cos  6 

- sin  <P  sin  \p 

+ cos  <{>  cos  \p  cos  0 

cos  \p  sin  0 


sin  $ sin  0 


- cos  <f>  sin  0 


cos  0 


We  note  that  each  element  of  SO (3)  can  be  parametrized  by  some  triple  of 
Euler  angles;  however,  when  0 * 0 or  it  this  parametrization  is  not 
unique. 

To  obtain  useful  special  functions  that  are  defined  on  SO (3)  we  shall 
consider  the  matrix  elements  of  a sequence  of  group  representations  of  S0(3). 
It  should  be  recalled  that  a group  representation,  which  is  described  in 
detail  in  [5]  - [7]  can  be  thought  of  as  simply  a group  of  matrices  to  which 
the  group  S0(3)  is  homomorphic.  In  particular,  we  will  use  the  sequence  of 
finite-dimensional  unitary  representations  {D  (4>,0,^)»  i “ 0,1,...}  attri- 
buted to  E.  P.  Wigner  whose  components  are  described  in  [6,  p.144]  by: 


Q .ra-n  -inx|>,J l /n.  -imp 

D - i e ^d  (0)e 
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(2) 


(3)  d£  (0) 
ran 


sinn~mQ (14cq3  9)m ( &-n)  1 ^ 

2l[  a+m)  ! a-®)  l l/™  *. 


,Z+n 
a 

d(coa 


(cos  0 - l)^'+m(l+cos  0)*'_m  , 


where  m and  n are  integers  such  that  - Z * m,n  i Z . 

o a Z 

It  is  proved  in  [7,  pp. 233-284]  that  the  functions  D~(R)mn  ■ D (<{>,0,^)^ 
form  a complete  orthogonal  system  in  the  space  of  square  integrable  functions 
f((j),9t’^)  , 0 £ 0 < tt  , 0 < Q,\l)  < 2tt  , with  respect  to  the  inner  product 


<fi’V 


i f2Tr  r r2ir 

Jo  J0 


fjC't’.e.*)  sin0d<})d0d^  . 


The  completeness  of  this  system  in  the  Hilbert  space  of  square  integrable 

f unctions  defined  on  the  rotation  group  is  analogous  to  that  of  the  circular 

functions  on  the  circle,  , and,  also,  that  of  the  spherical  harmonics  on  the 

2 

sphere,  S 

We  now  define  a rotational  exponential  Fourier  density  of  order  N to  be 
denoted  by  REFD(N)  as  a probability  density  defined  on  S0(3)  of  the  form 


(4) 


p (R)  = pOM.'W  =exP 


- i i 

£=0  m,n=-J6 


0 £ 
where  Sqq  is  a normalizing  constant  and  all  other  coefficients  a are 

arbitrary  complex  numbers.  In  addition  to  the  completeness  property, another 

reason  for  our  choosing  this  class  of  special  functions  is  that  for  continuous 

densities  we  have  the  following  approximation  theorem,  which  is  a generalization 

of  Theorem  1 in  [1], 
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Theorem  1.  Let  p be  a probability  density  on  SO (3).  Assume  that  p is 


continuous.  Then  for  any  given  positive  number  e , there  exists  an  REFD, 

N SL 

= e*P  I l » such  that  1 1 P “ PN I I < e • 

1=0  m,n=-£ 


Proof.  Assume  that 


inf{p(x)  : x e SO  (3)}  = c > 0 . 


This  assumption  can  be  removed  in  exactly  the  same  way  as  in  the  proof  of 
Theorem  1 of  [1],  We  note  that  f(x)  =SLn  p(x)  is  now  well  defined  and  also 
continuous  on  SO (3). 

Since  S0(3)  is  compact,  in  view  of  the  Peter-Weyl  Theorem  [6,  p.99],  for 


any  0 < 6 < 1 , there  is  a linear  conbination  of  D , say 
N(S)  l mn 

f<$  = I 1 amn(6)Dmn  » such  that  I I f6  “ f I I oo  < e • Xt  follows  that 

1=0  m,n=-£, 

llf6IL  < 1 + IlfIL  = sm  . 

Define  a function  g : R'*'  -*■  R’*"  by 


rexp  x , x £ M"! 

g(x)  = \ ^ 

l exp  M , x S MJ 


and  an  operator  G on  the  set  of  real  functions  on  S0(3)  by  (Gu) (x)  = g(u(x)) 

It  is  obvious  that  g satisfies  the  Caretheodory  conditions  [14,  p.20]  and 

2 2 

G transforms  every  function  in  L (SO (3))  into  a function  in  L (SO(3)). 

By  Theorem  2.1  of  [14,  p.22],  the  operator  G is  continuous.  Hence 
given  any  e > 0 , there  exists  a 6 > 0 such  that  if  | | f ^ - f| | < 6 , 

then  ||Gf^  - Gf| | < e . We  choose  f^  = f^  . Then  ||exp  f^-  p| | = 

| | Gf ^ — Gf | | < e . This  completes  the  proof. 
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Remark.  We  note  that  is  not  necessarily  the  N(6)-th  order  partial  sum 

SN(<5)^  oF  Fourier  expansion  of  f . Thus  in  contrast  to  Theorem  1 

of  [1],  it  is  not  necessarily  true  that  X.im||exp  s f - p|  | = 0 , which  is 

k-*» 

very  useful  in  determining  p^  in  the  statement  of  Theorem  1.  Such  a 

deficiency  can  be  removed,  if  we  require  p to  be  twice  continuously 

differentiable.  This  extra  requirement  ensures  that  Jin  p for  p > 0 , 

may  be  expanded  in  an  absolutely  and  uniformly  convergent  Fourier  series  in 

terms  of  the  rotational  harmonics,  [15,  p.513].  Thus  the  Peter-Weyl 

mn 

Theorem  is  not  needed  for  uniform  approximation  of  a continuous  function. 

A remark  is  in  order  about  (4).  From  our  choice  of  Euler  angle  para- 
metrization,  it  is  easily  seen  that  if  the  range  of  permissible  values  of 
<}>,0,ip  is  ignored,  then  a rotation  R'*'  with  Euler  angles  (<(>-HT,  - 0,(jH-Tr) 
would  be  equivalent  to  the  rotation  R with  Euler  angles  , 0 in  the 

sense  that  the  final  positions  of  the  coordinate  axes  would  be  identical. 
Consequently,  it  is  advantageous  to  extend  the  definition  of  the  function 
D ($,8,^)^  so  that  we  have  the  property 

(5)  p(R)  = p (R1) 


and  so  that  we  can  lift  our  restrictions  on  0,  <p,  and  ip,  except  0^0. 

£, 

First  we  extend  the  definition  of  d (0)  to  include  all  values  of  0 such 

mn 

that  0 < ( 0 1 < 7T  by  defining 

d£  (8)  = (-l)nH'nd£  (-0)  , - IT  < 9 < 0 . 


i/sing  (2)  we  see  that  (5)  holds.  For  the  situation  where 
shown  in  [6,  p.114]  that 


mn 


-im(<J)-Hj;)  r 

mn  ’ 


0 


0 , it  is 
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.**• ■ 


t. -mansm. 


thus,  since  <t>  + ijj  is  the  fixed  amount  of  rotation  about  the  z-axis  even 

though  each  angle  is  not  uniquely  determined,  ve  are  assured  that  for  any 

choice  of  Euler  angles  representing  this  rotation  as  well  as  any  other 

rotation  R the  density  (4)  is  well-defined. 

Before  proceeding  to  consider  an  estimation  problem  on  S0(3),  we 

£ £ 

enumerate  some  properties  of  the  functions  D (R)  and  d (0)  which  will 

mn  mn 

be  employed  in  Sections  IV  and  V: 


(6) 


D*«lVn»  ‘ l 

s=-£ 


sn 


for  any  rotations  and  R^  belonging  to  S0(3),  [7,  p.281]; 


(7) 


dl  (0)  - (-1)“^  (9)  = d*  (0)  , 16,  p.157] ; 

mn'  — m , — n — n , — m 


0 o.  . . A L\A  V L\  T 

(8)  D <R>D  00_ini  = eim4> l (2L+1>  ](  *>>« 

11111  m n L \m  m'  M/\n  n*  11/  ^ 


Imp 


where  M “ -m  - m'  and 
L for  which  i 

Wigner , 


N = -n  - n*  , [6,  p.157].  The  summation  is  over  all 
L £ H + l'  and  for  which  the  3-j  coefficients  of 


rl  L> 
tm  m'  11) 


and 


'l  V 


in  n1  N , 


* 


exist;  these  coefficients  are  defined  in  [6,  p.120]  by 


2.1-k+n  r q+k-A)  1 (k+SL-j  ) I U+i-U  8 (&+p)  ! (JE-p)  tl  */z 
L (j+k+fc+1)  ! (j+m)  ! (J-m)  I (k+n)  1 (k-n)  ! J 

r ( }t  (ft+J-n-t?  ! (k+n+t)  ! 

x 4 *■>  (Jl+p-t)  ! (t+k-j-p)  Itl  (£-k+j-t)  1 
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f i 


a.  -tv**.**. 


where  the  summation  is  over  all  integers  such  that  the  arguments  of  all 
factorials  are  nonnegative. 

III.  An  Error  Criterion  and  Optimal  Estimates  on  SO (3). 

In  order  to  define  an  error  criterion  for  orientation  estimation,  it 
is  necessary  to  have  a measure  of  the  distance  between  two  orientations. 

We  will  first  describe  such  a measure,  using  quaternions  [11].  We  recall 
that  a rotation  about  an  axis  in  the  direction  of  a unit  vector  [£,m,n]' 
through  an  angle  $ is  represented  by  the  (unit)  quaternion 

q = ^ql,q2,q3’<14],*=  ^C°S  2*  113111  2*  mSln  f’  nsln  fJ  ' 

Given  two  orientations,  the  minimal  angle  in  radians  required  to  bring 

one  into  the  other  is  a natural  measure  of  distance  between  them  and  defines 

a Riemannian  metric  on  S0(3).  If  the  orientations  are  represented  by  the 

quaternions,  q and  p , and  the  minimal  angle  is  denoted  by  p(q,p)  , 

then  we  have  q'p  = cos  |p(q,p)  . As  y(l  - cos  p)  is  a monotone  increasing 

function  of  p , a measure  of  distance  between  p and  q can  be  defined 
1 2 

to  be  -^(1  - cos  p (q,p))  = 1 - (q'p)  • It  can  be  shown  that  if  the 

orientations,  q and  p , are  described  by  the  3x3-dimensional  orthogonal 
matrices,  Q and  P , then  this  measure  of  distance  can  also  be  expressed 
as  ^-(3  - tr  PQ')  . 

We  are  now  ready  to  define  the  error  criterion  for  orientation  estimation. 
Let  q be  a random  quaternion  and  p its  estimate.  Then  a measure  of  the 
estimation  error  is 

(9)  J(q.p)  = E(1  - (q'p)2)  . 
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If  the  probability  distribution  of  q is  given,  the  estimate  p which 


minimizes  J may  be  obtained  from  observing  that 


J(q,p)  = 1 - p'  E(qq')p  . 


It  is  well-known  [12,  p.142]  that  the  quadratic  form  p’Vp  of  the  positive 
definite  matrix  V = E(qq')  is  maximized  when  p is  the  unit  eigenvector 
associated  with  the  largest  eigenvalue  X of  V . Moreover,  the  maximum  value 
is  X . Hence, 

min  J(q,p)  = 1 - §'E(qq')§ 

P 

= 1 - X 

where  X = the  maximum  eigenvalue  of  E(qq') 

q = the  unit  eigenvector  of  E(qq') 
associated  with  X . 

The  probability  distributions  on  SO (3)  are  expressed  in  terms  of  the 
Euler  angles  (cf> , 0 ) in  this  paper.  The  following  relationships  [10,  p.380] 
between  the  quaternion  components  and  the  Euler  angles  will  have  to  be  used 
to  calculate  the  optimal  estimate  q and  its  estimation  error  1 - X : 


(10) 


0 ±hb 
q^  = cos  — cos 

q2  = sin  — cos 


q3  = ’in  2 sin  V 


0 . 

q^  = cos  — srn 


Once  vector  q is  determined  we  can  immediately  determine  the  Euler  angles 
(5> , 0 , ip)  for  the  optimal  orientation  estimate  from  the  relations: 
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(11) 


COS0  = 2(5l2  + q42)  - 1 , 
sin$  = - j (q.^  + q2q4)  , 


sin$  = j (q2§4  - q.^)  , 


A = / 2 a.  2.  .a.  2 a 2. 

(qi  +q4  )(q2  +q3  > 


A 

0 * 0 i IT 


cos$  = J i^2  - §3q4) 

C03$  = J (^^2  + §334) 


which  are  the  inverse  of  the  above  set  of  relationships  (10). 


IV . Estimation  for  Rotational  Processes  with  Rotational  Noise. 

Suppose  s is  a rotation  of  a rigid  body  with  Euler  angles  (<J>,0,4') 
and  v is  a second  rotation  with  Euler  angles  (e,ri,6)  . We  can  obtain 
a product  rotation  m with  Euler  angles  (<J>,6,^)  by  the  successive  rotation 
of  the  rigid  body  by  s and  then  v . We  denote  this  product  by 

(12)  mOjj.Q.'f')  = v(e,n,S)  • e(<p,9tip)  . 

Our  reason  for  writing  v to  the  left  of  s is  to  be  consistent  with  the 
matrix  equation  resulting  by  the  use  of  the  orthogonal  matrix  representation 
where  the  matrix  S representing  s is  pre-multiplied  by  the  matrix  V 
representing  v to  obtain  M = VS. 

Performing  the  indicated  multiplication  when  each  matrix  has  the  form  (1) 
and  equating  respective  elements,  we  may  obtain  nine  relationships  among  the 
Euler  angles  of  m,  v,  and  s,  which  uniquely  determine  (<p, 9,ip)  in  terms  of 
(e  ,n,5)  and  (<p,9,ip)  . The  relationships  are  very  cumbersome  and  thus  will 
not  be  displayed  here. 

We  now  consider  a rotational  signal  process  which  is  a sequence  of  rotations 
(sk}”_0  defined  on  S0(3)  whose  terms  are  related  by  the  equation 
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(13) 


w,.  ^ 


‘Hc.+l  "k' 

where  is  a deterministic  rotation  whose  Euler  angles  are  known.  Let 

im.  ) and  {v,  } be  sequences  of  rotation  on  S0(3)  that  constitute  a 
k-1  k-1 

measurement  process  and  a measurement  noise  process,  respectively,  such  that 


(14) 


“k 


vk  0 sk 


The  first  estimation  problem  to  be  solved  can  now  be  stated  as  being  that 
of  finding  the  optimal  estimate  of  s^  given  the  set  of  measurements 


k A 


m 


- {m^,...,m^}  , k 1,2, 


using  the  error  criterion  (9).  From  our  results  in  Section  III  we  recognize 
the  fact  that  in  order  to  determine  the  optimal  estimate  which  involves 
conditional  expectations,  it  is  crucial  to  be  able  to  produce  the  conditional 
density  p(s^|m'C)  , which  we  note  can  be  calculated,  by  Bayes  rule  from 


(15) 


p(s  Im*)  =■  Ckp(m1(.|sk)p(sk|mk_1) 


where  Ck  is  a normalizing  constant. 

It  will  now  be  demonstrated  that  REFD's  introduced  in  Section  II  are 
ideal  ones  to  use  in  computing  this  conditional  density.  Suppose  Sq  and 
{ vk)  have  REFD(N)  (if  they  have  different  orders,  we  can  let  N be  the  maximum 
order  and  by  inserting  zero  coefficients  make  all  densities  of  order  N)  which 
are  described,  respectively,  by 


(16) 


(17) 


N l „0  , 

p(sQ)  - exp  I l D (sQ) 

£■0  m,n=-x, 


inn 


N H „ 

p(v,)  - exp  l l b*k  D^V.) 

K on  L 0 “n  k 'mn 

Jl-0  m,n— £ 
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Ic— 1 

Let  us  further  assume  that  the  conditional  density  pCs^Jm  ) is  a 


REFD(N),  denoted  by 


N £ 


/ I k-ls  v V J*-/c  t 

(18)  PCSkrl|  in  ) = «p  I l a ° ^k-^mn  * 

X/**u  in  9 n=-x, 

k 

Under  these  assumptions  we  will  show  that  p(s^|  m ) is  also  a REFD(N) 

and  at  the  same  time  exhibit  a recursive  formula  for  the  Fourier  coefficients 


from  (13)  and  (14)  and  the  group  property  of  S0(3),  vfc  and  s^_^  can 


be  expressed  as 


vk  * “k  0 Sk1  and  sk-i  = wk-i  ° 8k 


so  that,  using  (18),  p(SjJm  ) is  a REFD: 


P(sk|mk  X)  = I l 1 ^"k-i  ° sk>, 

1= 0 m,n=-£ 


V v £,k-l  r n£,  -1  x J*,  x 

= exp  ) ) a ) D (w  , ) D (s,) 

v nLn  L „ mn  u 0 k-1  ms  k sn 

£-0  m,n=-£  s=-£ 

- - jo  ,.L  U ^ dS'-i)J  D‘(9t>“ ! 


for  the  second  equality,  property  (6)  was  used  while  the  third  equality  is 


obtained  by  a change  in  summation  indices. 

Likewise,  the  conditional  density  p (m^ | sfc)  is  a 


REFD: 


pOnjJV  - ~P  L * b*k  DA(mk  o s^1), 

£*0  m,n=-{. 


-jo  [ } , bj"  J Dt(^1>-' 

£=0  m,n=-£  uj  — £ J _ 

{ l [<-»"“  i ‘5%  “*<%>!,-] "‘“k'-,  • 

-0  m,n«-£  L j— * J*  J 
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1 


the  last  equality  results  from  a change  of  summation  Indices  and  the  fact  that 
if  has  Euler  angles  (4>,0,^)  then  s^  has  Euler  angles  (tthJj,0,it-<J>) 

so  that,  using  (2)  and  (7)  , 

/ -lx  ,m-n,  . xin+n  lm iKf.  x in6 


dV"1)  - (0)ein4) 

k ran  mn 


, , x nrHi . -n+m  imiK&  ,Qx  ind> 
= (-1)  x e rd  C0)e  Y 
— n,  — ra 

- (-l)®fnD*'(s  ) 

k -n,-m 

Substituting  Cl9)  and  (20)  into  (15),  we  obtain  the  REFD 


p(sk|m  ) = Ck  exp 


-f  z [! 

£-0  m,n=-£  Lj—  £ 1 1 J1 

_l)mtn  fe£k  D£  ( j \1  D*  (s  ) 

j,-m  k j,-n J J kmn 


Thus  p(s^|m^)  is  a REFD(N)  whenever  p(s^_^|m^  ^)  is  a REFD(N) 


but  since 


P(s0IV  " P<s0> 


is  itself  a REFD(N)  by  induction  we  have  proved  the  following: 

Theorem  2.  If  {s^}  is  a sequence  of  random  rotations  on  SO (3)  where 

si  “ W.  o s, 

k+1  k k 

for  some  sequence  of  deterministic  rotations  and  m^  is  the  corresponding 
measurement  sequence  where  each  sfc  has  been  corrupted  by  a noise  rotation  vk 


such  that 


"k  = vk  0 sk 
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£dfc>  * i 


"t  \f..-  . 


then  if  Sq  and  {v^}  have  REFD's  given  by  (16)  and  (1?) > respectively, 
then,  for  k = 1,2,...,  the  conditional  density  p(s^|m  ) is  a REFD(N) 
of  the  form 


p(sk|mk)  - exp  I l D*(s  ) 

£=Q  m,n=-£ 


mn 


ilk. 

where  the  coefficients  a^  are  determined  recursively  by  the  formulas 


ilk 

a 

mn 


i— £ 


L jn  v k-1  jm  j,-m 


£ t 0,  k = 1,2,3, ... ; 


and  a^Q  is  a normalizing  constant. 

This  theorem  enables  us  to  calculate  the  sequential  conditional  densities 
by  updating  recursively  a finite  and  fixed  number  of  parameters.  Using  the 
conditional  density  p(SjJmS  at  time  k , the  optimal  estimate  of  the 
orientation  can  be  determined  as  suggested  at  the  end  of  Section  III.  Namely 
we  first  compute  the  conditional  covariance  matrix  E(q (k)q'  (k).|*“k)  where 
q (k)  is  the  quaternion  for  s^  whose  components  expressed  in  terms  of  the 
Euler  angles  are  given  by  (10) . Then  we  use  some  standard  numerical  method 
such  as  the  power  method  [13,  pp. 147-150]  to  compute  the  largest  eigenvalue 
A(k)  and  the  associated  unit  eigenvector  q(k|k)  . The  Euler  angles  of  the 
optimal  estimate  may  then  be  determined  from  q(k|k)  through  (11).  Tire 
estimation  error  is  1 - A(k)  . 


V.  Estimation  for  Rotational  Processes  with  Additive  White  Gaussian  Noise. 

A second  model  for  v/hich  the  estimation  problem  can  be  solved  using  REFD's 
is  described  by  the  equations 


(21) 

8k+l  " wk  * 8k 

(22) 

mk  - h(sk)  + vk 
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’/here  {a^}  Is  again  a sequence  of  random  rotations  belonging  to  S0(3)  and 
{v^}  is  now  a sequence  of  p-dimensional  vectors  of  observed  outputs  of  the 
signal  process,  h is  a vector-valued  function  defined  on  SO (3)  that  is 


is  a sequence  of  p-dimensional  independent 


Gaussian  vectors,  each  with  zero  mean  value  and  covariance  matrix 


The  completeness  property  of  the  functions 


any  e > 0 and  for  each  component  hJ  of  function  h 


Each  noise  vector  has  density 


then 


% ■ W + \ 


(23) 


'<MV  ■ «’')"p/2(“«  %)'1/!^p  -ft  f ■£*  k - ? I 

v i j j “1  L X(®0  in  ^ iv"  “ X*  m 

< R - f I h£j  D£(0  1} 

l>  £-0  m,^-H  13111  kmnJJ 

- (2k)-»/2Cd«t  v‘y'«p  (c0  + X ! tclAv„, 

^ Jc»0  m,n®-x 


M l V not  o 

+ I I I c (m,n,m'  ,n'  )D  (s.) 

£,£'=()  m,n=-H  m',n 


where 


co  = -i  ^4j  mW 


cl  = ~ f R^lml  h£j  + mi  h*1] 

mn  2 . . , ic  |_  k mn  H mnj 

JU',  , 1 ? ij.Hl  .H'j 

C (m,n,m  n ) = - j l R,Jh  h ,\ 
z . . - k mn  m n 


But  properties  (7)  and  (8)  imply  that 


<“>  (n  m'  iXn  u'  »)  DV-«,-» 

k-1 

so  that  the  expresion  given  in  (23)  is  a REFD(2M).  Now,  if  pCs^  jj  m ) 
is  a REFDCN^)  then  by  (15)  p ( m^)  is  a REFD(max{2M,N^  ^})  , so  if 
has  a REFD(N)  given  by  (16)  then  p(SjJm  ) is  a REFD(Max{2M,N})  for  any  k , 
with  a finite  number  of  recursive  formulae  existing  for  the  coefficients. 
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These  formulas  will  not  be  enumerated  here  since  they  are  tedious  to  list  for 


an  arbitrary  integer  M . 

Ue  exemplify  this  procedure  for  p * 3 where  the  vector-function 
h has  a very  simple  form  so  that  the  computations  are  not  unduly  .Cumbersome: 


h(sk)  = 


cos  9, 
k 

cos  <f>k 

COS 


sin  0. 
k 

sin  9. 
k 


the  components  being  three  of  the  direction  cosines  of  the  matrix  representation 
of  sk  of  the  form  (1). 

Let  us  assume  that  {vk}  is  a vector  Gaussian  process  with  each  term 
v^  having  zero  mean  and  3x3  covariance  matrix  Rk  , and  that  Sq  has  the 
REFD(N)  (16).  Now 

pCmjJ  sIc)  = (det  R^  ^exp  - jfm^.  - h(sk)]'  R^m^  - Ms^)]  . 


Using  (2),  (3),  and  (8)  in  the  following  calculations,  we  first  write  h(sk) 
in  the  form 


h(sk)  = 


»X> 

00 

i[D1(sv)  + D1  (s.  ) ] 


Si 


10 


-1,0 


-!D1(sk)  + D1(s  ) ] 

Si  oi  k 0,-1 


so  that 


_ 19  _ 


■h(\),Vvh(V)  • 

+ d1(V-i,o)1^2+2<Vd1(VooiiVi|(d1(Voi+d1<Vo,-iK3 
+ Iv|(I>1(8k)io+I,1(8k>-i,o)l2^2+2[vJ<cl(sk)io+Dl<9k)-i,o)l 

X lvi<D1(sk)01«1<sk)0i_1]l^3+[^-i(D1(kk)01+D1(8k>0,-l)l2^3 


■ j.  .^Vk  - 2 j.^Voo-  n 

i,j=l  j-1  j-1 

- •K  ^ 1 j ■i'kV<8k>o,-i 

j-1  j-1  j-1 

+ I Ri1|D°<8k)00+2I)2(8k)001t§  1 

+{l 1 ^‘^kV-XVoil-N;2'  J °2<8k>20+  7">2(8k)00 

- D°<8k)00,+>f  D2(8k>-2,0>-  J ^II>2(8k)U-Dl(8k>U+Dl(8k>l,-l 

+I>2(8k)i,-i+1>1(8k)-i,i+D2(8k>-i.x+I,2<8k)-i>-rDl(8k>-i,-i1 
' \ I/=°  <8k)02+ I <B  *VoO~D  <8k)00)+/J  0 *8k30,-2^‘ 


Assuming  that 


p (ek— 1 1 n*k_ 


v r k-1  1 

exp  l l a * D (a.,) 


1-0  m,n--l 


k-1  mn 


? 


and  using  (15)  and  (19)  we  obtain 


N i 


p(s  jm*)  = exp  I l DA(sk) 

1 1=0  m,n— Jt 


Ok 

where  a^^  is  a normalizing  constant. 


lk  o v j V l.k-l-l,  -1  x 

aoo  " ' 2l1  \ \ + jL/JO  ° (,WjO  * 


? J V -1 


j-1  k k i=-l  jn 


k-1  jm 


lk  1 „23 


l . V l,k-l_l,  -lx  ll  , 

+ ,L1aj»  “ (’Vi)jm  , M ■ 1 


2k  _ l,ODll  „22  33N  r x ,.R-x  x 

a00  3(2Rk  _ Rk”  Rk  * ,_S  jO  ° ( k-l^jO 


2,. k-1  2.-1 


2k  1 r-PP  , v 2,k-l_.2/  -lx  n |_.  | 0 

amn  = ^ + lajn  D (wk-l>lm  , «®  = 0 , | mta | = 2 , 


p = 2 


5(!»|,2)x,5(|n|,2)  . 


fi  . _lp  2,k-r  r 2, k-1  2,  -1  N - 

!V3  1 \ + am  + .L2ajn  ° ^k-Pjm  » = 0 » 


1 , p - 25(|m|,l)x36(|n|,l)  . 


2k  1 „33 


- 1 + JL2aj^1,,2<vi)j-  * IH  ° 1 * |m|  + |n|  “ 2 


«,k  r J.k  JL,  -lx 
a = ) a.  D (w.  ,).  , otherwise, 

ran  jn  k^l'jm 


[We  have  used  the  symbol  6(a,b)  to  denote  the  number  1 when  a = b and  0 
otherwise].  Now  optimal  estimates  can  be  obtained  exactly  as  in  the  last  section. 
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VI.  Detection  for  Rotational  Processes 


In  this  section  we  consider  the  detection  problem  of  determining  the 
presence  or  absence  of  a rotational  signal  process  using  each  of  the  measure- 
ment processes  given  in  the  past  two  sections.  It  will  be  shown  that  the 
likelihood  ratio  of  the  presence  to  the  absence  of  the  signal  can  be 
efficiently  ascertained  when  REFD's  are  utilized. 

Let  us  first  refer  to  the  rotational  signal  and  measurement  processes 
{s^}  and  {m^l  described  by  (13)  and  (14)  where  s^  and  (v^l  have  REFD's 
given  by  (16)  and  (17)  which  are  independent  processes.  We  now  consider  the 
following  hypotheses  that  the  signal  rotation  s^  is  present  and  absent 
respectively: 

**i  ; “k  ■ \ • sk 

Ho  ! \ ■ \ 

The  likelihood  ratio  that  we  want  to  compute  is  given  by 

,,  k.  4 P<»k!Hl> 

L(“  } = kl 

P(®  I H0) 

Since  (s^)  and  {v^}  are  independent,  we  have 

pOn6!^)  = E[p(mk|  H1,  sk)(  Hj 

k « l £j  _ 

= E[  II  exp  y y b D (m  o s ) ] 

1 v L L mn  j j mn 


j»l  £,=0  m,n=-£ 


while 


N l 


p(n^|H  ) = n exp  J l b^  Am ) 
j=l  £=0  m,n— 1 J 
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k*  * * - 


so  that 


where 

(26) 


sT1)  - 0*001 

j mn  j mn 


k 

= E[exp  l W(s  ,j ,m  ) ] 
j = l 3 3 


W(s  ,j,m  ) = l l b^[D£(m  . s"1) 

2 3 4- 0 m,n— £ 3 3 


- D (m  ) ] 

mn  j mn 


We  now  state  a lemma  whose  proof  is  identical  to  that  of  a corresponding  lemma 

proved  in  [1,  pp. 14-15]  with  the  circular  function  replaced  with  the  functions 

D*(s)  : 

mn 

Lemma  1.  Let  { s^}  be  a random  signal  process  on  S0(3)  and  {v^}  a white 
noise  process  on  SO (3)  which  is  independent  of  { s^}  , having  the  rotational 
exponential  Fourier  density  (23).  Let 


T * { t ^ , • • « , t^} 

T'  = { t , • . . ,t^,} 

Sq,  = {st  » • • • >s^  } 

i q 

where  q,q',t^,  and  t|  are  positive  integers.  If  the  measurement  process  is 
described  by  mk  = vk  0 then  the  conditional  joint  probability  density  is 

k E[exp  H(l,k) lSTtST, ]p(ST  |ST* ) 

p(Sr|m  ,ST,)  - E[exp  H(l,k) |ST,] 
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r. 


where  R(l,k)  = £ W(s  ,j,m  ) and  W(s.  ,j,m  ) i 

j=^  J J J J 


is  given  by  (26) . 


We  can  now  obtain  a recursive  formula  for  L(m  ) by  using  the  smoothing 
property  of  a conditional  expectation  and  Lemma  1: 


k 

L(m  ) = 


IV 

E[exp  £ W(s^,j,mj)|Sk  ]p(sk)dsk 


k-1 

^ E [ exp  l W(s^,j,nu)jsk]exp  W(sk,k-,mk)p(sk)dsk 


k-1 


k.-l. 


E[exp  £ W(s..,j,m.)]p(sk|m  ' )exp  W(sk,k,mk)dsk 

j=l  3 3 


E[exp  H(l,lc-1)] 


p(sk|mk  1exp  WCs^.k.ir^Jds^^ 


= L(mk  ^QOnS 


C . r2ir  tt  r 2tt 

where  J f (sk )dsk  = J ^ »\^sin  ^k^k^k^l 


and 


k.  A 

Q(m  ) = 


k-1 


P (sk I m )exp  W(sk,k,n^)dsk 


The  integral  Q(in  ) can  be  evaluated  in  the  following  form  which  is 
obtained  with  the  assistance  of  (19),  (26)  and  (6): 

N l 


0 (nk  ) 


exp[W(sk,k,mk)  + l l a^k_1  D^w^o 
£=0  m,n=-L 

r r r , L,k-1  -1  v nSL,  v 

exp  I I I {a  D (vk  ) D (s  ) 

L 1-0  m,n=-£  j=-£  J J 


+ b>VvV>3„  - D"'<Vmn1)dsk 
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r warn* 


N £ £k  5 

exp{  “ I 1 b D (m.  ) } 

r nLn  L „ mn  k nrn 
£=0  ra,n=-£ 


<p  j?  £ C£k  D£(sv)  ds, 

. L . L . mn  K mn  k 


j £=Q  m,n=-£ 


where 


C^k  = l {aJk  D^w"1  ) + (-l)m+tl  bf"  Atm).  } 

mn  ^_£  jn  k-l'jn  j,-m  Tt  j,-n 


Let  us  now  study  the  detection  problem  that  results  from  using  the 
signal  process  discussed  in  Section  V described  by  (21)  and  (22). 
Suppose  we  have  the  hypotheses 

Hi  : mk  = h(sk)  + vk 


H0  : mk  = \ 


Using  (23)  it  is  found  that  the  likelihood  ratio  is 


L(m  ) = E[exp  H(l,k)] 


11(1, k)  = l W(sj,j,mj) 


j=l 


M £ ,£  . v v v'  JM' 


W(s  j,m  ) = I l C D Cs  ) + . I,  . I _ . 1.  _c  (m,a,rn',n’) 

J J £=0  m,n=-£  J 


£,£'=0  m,n=0  m',n,=0 


x As ,)  D£'(s  ) , , 
j mn  3 m n 


„£  1 ? „ijr  j i&j  * i ,£i, 

C = ■=•  ) R.  J [m,  h J + m,  h ] 

mn  2 " , k k mn  l mn 


rW , , ,»  _ 1 v .«■  .*'j 

C (m,n,m  ,n  )■*-■=-  / R,  h h,, 

* ’ ’ 2 (•_  k mn  m n 

1 > J J- 


An  argument  identical  to  that  used  in  the  previous  detection  problem 
yields 
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F 


^ ■ r 


L(pik) 


El  exp  H(l,k)|  sk]p(sk)dsk 
L(jnk  X)  p(sk|mk  1)exp  W(sk>k,rak)dsk 


= L(mk  1)Q(mk) 


Since  it  has  been  shown  in  the  last  section  that  the  conditional  density 
k-1 

p(sk|m  ) is  a rotational  exponential  Fourier  density  of  the  form 

N £ 


v r £k  _£,  . 

exp  l l a D (s.  ) 
o a o ran  k ran 

>c=»0  m,n=-£ 


where  the  coefficients  can  be  computed  recursively,  we  can  write  the  integral 
k 

Q(m  ) in  the  form 


. M £ 0 0 

Q(j^)  = exp  £ l C D (s  ) 
JZ  £=0  m,n =-£  n k 


v v ZZ'  Z £' 

+ I I l C (m,n,m*  ,n'  )D  (s  ) D (sfc ) , , 

(j  n f ntin  k ran  k m n 

X/,36  =0  m,n=-Jc  m ,n  =-£ 

» * £k  l 

+ l l a D (s.  ) d3  . 
oLn  L a mn  k mn  k 

x.=0  m,n=-£ 


Consequently,  Q(m  ) is  of  the  form 


P £ 

exP  1 l a?k  D*fe  ) 

\ £=0  ra,n=-£  mn  ^ k. 


ds, 
mn  k 


£k  £k 

where  P = max{N,2’t}  and  the  coefficients  a which  are  functions  of  a 

mn  mn 

£ £i 

^mn  * ^mn  ’ ’ anc*  wk  are  °^ta^-ne<^  using  (24).  Hence  the  computational 

scheme  for  the  likelihood  ratio  L(m  ) is  finite  dimensional.  Rather  than 
attempting  to  exhibit  the  formulas  in  the  general  case  we  will  produce  them  for 
the  simple  example  at  the  end  of  Section  V where 
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hO*) 


cos  <t>k  sin  0^ 
cos  sin  0k 


and  (v.  } is  the  Gaussian  with  zero  mean  and  covariance  matrix  • 

can  be 


Since  it  was  shown  that  the  conditional  density  p 
written  as 


(sj  m^) 


N 4 


4k  „4, 


p(s Jm1^)  * exp  l l D (sk\nn 

K 4*0  m,n=-4 


with  the 
ratio  is 


coefficients  given  by  (25),  it  is  easily  verified  that  the  likelihood 


L(mk)  = E[exp  H(l,k)] 


H(l,k)  = l W(s  ,j,m  ) 

j=l  3 


N 4 


w(a  .j.-j)  - i i & - K IS  ■’  "5 

3 3 £_o  m,n=-4  s,t-i 


and 


(mS  = exp{a°^+  a°’k  1-  l I \ % 


l00  00 
N 4 


s,t=l 


exp 

^ 4*1  ra,n=-4 


l l .Cmn  D*'(sk)mn“"k 


ds. 


C^-a^t  l a^D^w  ) 
mn  mn  jn  K 1 3 
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VII.  Conclusions 


In  this  paper  we  have  formulated  and  solved  estimation  and  detection 
problems  for  discrete-time  processes  that  arise  on  S0(3).  The  signal 
process  s^+^  is  assumed  to  be  obtained  from  s^  by  a concatenation 
which  is  the  successive  rotations  s^  and  w^  where  the  latter  is  a known 
rotation,  while  two  types  of  measurement  processes  have  been  used:  (i) 
the  observation  m^  is  the  concatenation  of  and  a noise  rotation  , 

and  (ii)  the  observation  m^  is  a smooth  function  of  s^  corrupted  by 
additive  white  noise. 

An  error  criterion  which  differs  from  a least  squares  criterion  such 
as  appears  in  [8]and  [9]  has  been  presented  together  with  the  resulting 
optimal  estimates.  In  addition  a new  density  function,  the  rotational 
exponential  Fourier  density  has  been  introduced  which  can  be  used  to  approxi- 
mate probability  densities  on  SO (3)  that  are  continuous. 

Since  this  class  of  density  functions  is  closed  under  the  operation  of 
taking  conditional  densities,  we  have  been  able  to  obtain  recursive  schemes 
for  optimal  estimation  and  detection  for  a rather  large  class  of  problems. 
However  recursive  schemes  for  the  analogous  problem  when  the  driving  term 
w^  is  a stochastic  process  have  not  been  resolved  since  these  densities  do 
not  have  the  property  of  being  closed  under  convolution.  Perhaps  a different 
representation  for  densities  on  S0(3)  or  a different  model  of  the  signal  process 
can  be  found  that  will  yield  a recursive  solution  to  this  problem. 

It  is  believed  that  the  procedures  of  this  paper  can  be  extended  to  SO(n) 
for  n > 3 by  using  the  appropriate  special  functions  defined  on  SO(n)  that 

H 

are  analogous  to  the  functions  D on  S0(3). 

mn 
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